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Abstract. The data analysis of current Cosmic 
Microwave Background (CMB) experiments like 
BOOMERanG or MAXIMA poses severe challenges 
which already stretch the limits of current (super-) 
computer capabilities, if brute force methods are used. In 
this paper we present a practical solution to the optimal 
map making problem which can be used directly for 
next generation CMB experiments like ARCHEOPS and 
TopHat, and can probably be extended relatively easily 
to the full PLANCK case. This solution is based on 
an iterative multi-grid Jacobi algorithm which is both 
fast and memory sparing. Indeed, if there are Aftod data 
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work we are presenting focus on the first of these issues, 
namely the map-making step. 

Up to the most recent experiments, maps could be 
made by a brute force approach amounting to solve di- 
rectly a large linear problem by direct matrix inversion. 
Nevertheless the size of the problem, and the required 
computing power, grows as a power law of the data set 
size, and the limits of this kind of method have now been 



points I along the one dimensional timeline to analyse, 
the number ot operations is ot 0[J\'t d lnyVt 0( j) and the 
memory requirement is 0{Ntod)- Timing and accuracy 



reached (Borrill 2000). Whereas the most efficient develop- 
ment in this massive computing approach, i.e. the MAD- 
CAP package (Borrill 1999) has been applied to the recent 
BOOMERanG and MAXIMA experiments (|de Bernardis 



et al. 2000; Hanany et al. 2000) some faster and less con 



suming solutions based on iterative inversion algorithms 



issues have been analysed on simulated ArtCrtrijOir'S and 
TopHat data, and we discuss as well the issue of the joint 
evaluation of the signal and noise statistical properties. 



have now been developed and applied too (Prunet et al 



2000). This is in the same spirit as (Wright et al. 1996) 
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1. Introduction 

As cosmology enters the era of "precision", it enters si- 
multaneously the era of massive data sets. This has in 
turn showed the need for new data processing algorithm. 
Present and future CMB experimen ts in particular f ace 
some new and interesting challenges ( Bond et al. 1999 ) . If 



we accept the now classical point of view of a four steps 
data analysis pipeline : i/ from time-ordered data (TOD) 
to maps of the sky at a given frequency, ii/ from frequency 
maps to (among others) a temperature map, iii) from a 
temperature map to its power spectrum Ci, iv/ from the 
power spectrum to cosmological parameters and charac- 
teristics of the primordial spectrum of fluctuation, the ul- 
timate quantities to be measured in a given model. The 



and we definitely follow this latter approach. 

Design goals are to use memory sparingly by handling 
only columns or rows instead of full matrixes and to in- 
crease speed by minimising the number of iterations re- 
quired to reach the sought convergence of the iterative 
scheme. We fulfill these goals by an iterative multi-grid 
Jacobi algorithm. As recalled below, an optimal method 
involves using the noise power spectrum. We have thus 
investigated the possibility of a joint noise (and signal) 
evaluation using this algorithm. 

Section || presents in detail the method and its imple- 
mentation, while section |^ demonstrates its capabilities 
by using simulations of two on-going balloon CMB exper- 
iments, ARCHEOPSH flBcnoit et al. 2000| ) and TopHat^ 
The results are discussed in section |4|, together with the 
problem of the evaluation of the noise properties, as well 
as possible extensions. 



2. The method 
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2.1. Optimal map making 



The relation between the sky map we seek and the ob 

s erved data stream may be cast as a linear algebra system _ [^ T ]V _1 A] 1 A T N~ 1 



the least square solution which was used to analyse the 
"COBE" data ( [Janscn and Gulkis 1992| ), 



( [Wright ct aL 1996| [Tegmark 1997|) . Let t and p indices de- 



(8) 



note quantities in the temporal and spatial domains, and 
group as a data vector, dt, and a noise vector the tempo- 
ral stream of collected data and the detector noise stream, 
both of dimension Aftod- We then have 



A tp x p + n t , 



(1) 



where A tp x p is the signal vector given by the observation 
of the unknown pixoliood cky map, xp, which hao boon ar 
ranged 



In this paper we will actually deal only with this estima- 
tor. Nevertheless as a next iteration in the analysis pro- 
cess, we could incorporate various theoretical priors by 
expliciting V{x). For example, it is often assumed a Gaus- 
sian prior for the theory, i.e. V(x) cx exp(—XpC pp ,x P '/2) 
where C pp > = (x p Xp,) is the signal 
In that case the particular solution turns out 



p pp' 
covariance matrix. 

to be the 



ao a vector of dimonoion Nm±. The Af t 



Mr, 



"obser vation" matrix A therefore oncompaoooo the acan 
ning strategy and the beam pattern of the detector. 

In the following, we restrict to the case when the beam 
pattern is symmetrical. We can therefore take x p to be a 
map of the sky which has already been convolved with 
the beam pattern, and A only describes how this sky is 
being scanned. For the total power measurement (i.e. non- 
differential measurement) we are interested in here, the 
observation matrix A then has a single non-zero element 
per row, which can be set to one if d and x are expressed 
in the same units. The model of the measurement is then 
quite transparent: each temporal datum is the sum of the 
pixel value to which the detector is pointing plus the de- 
tector noise. 

The map-making step then amounts to best solve for x 
given d (and some properties of the noise). We shall seek 
a linear estimator of x r , 



Gispert 1998): 



Wiener filtering solution (Zaroubi et al. 199q ; Bouchet and 
Gispert 1996 ; Tegmark and Efstathiou 1996|; Bouchet and 



W = [C- 1 + A T N- 1 A] _1 A T N- 1 . 



(9) 



<p i 



W pt d t 



(2) 



To motivate a particular choice of the Af P i X x Aftod matrix 
W , a Bayesian approach is convenient. Indeed we are seek- 
ing the optimal solution to this inversion problem which 
maximises the probability of a deduced set of theory pa- 
rameters (here the map x p ) given our data (dt) by max- 
imising V(x\d). Bayes' theorem simply states that 



V(x\d) 



T{d\x)T{x) 
V{d) 



(3) 



If we do not assume any theoretical prior, then x follows 
a uniform distribution as well as d. Therefore, 



But this solution may always be obtained by a further 
(Wiener) filtering of the COBE solution, and we do not 
consider it further. 

The prior-less solution demonstrates that as long as 
the (Gaussian) instrumental noise is not white, a sim- 
ple averaging (co-addition) of all the data points corre- 
sponding to a given sky pixel is not optimal. If the noise 
exhibits some temporal correlations, as induced for in- 
stance by a low- frequency 1// behavior of the noise spec- 
trum which prevails in most CMB experiments, one has 
to take into account the full time correlation structure of 
the noise. Additionally this expression demonstrates that 
even if the noise has a simple time structure, the scan- 
ning strategy generically induces a non-trivial correlation 
matrix [A T A^ _1 A] of the noise map. 

Even if the problem is well posed formally, a quick look 
at the orders of magnitude shows that the actual finding of 
a solution is non trivial task. Indeed a brute force method 
aiming at inverting the full matrix [yl^iV -1 ^!] , an op- 
eration scaling as 0(Af pix ), is already hardly tractable 
for present long duration balloon flights as MAXIMA, 
BOOMERanG, ARCHEOPS or TopHat where Af to d ~ 
10 6 and Afpix ~ 10 5 . It appears totally impractical for 
PLANCK since for a single detector (amid 10s) Aftod ~ 10 9 
and Afpi X ~ 10 7 ! 

One possibility may be to take advantage of specific 
scanning strategies, and actually solve the inverse of the 



convolution problem as detailed in (Wandelt and Gorsk 



T(x\d)\ rr P(d\r) 



4^- 



If we further assume that the noise is Gaussian, we can 
write 



2000). This amounts to deduce the map coefficients ai rn 
in the spherical harmonic basis through a rotation of a 
Fourier decomposition of the observed data. The map will 
then be a simple visualisation device, while the ai m would 



cx exp(—n t N tt , nc/2) 


(5) 


in ( 


3ouchct and Gispert 1996; Tegmark and Efstathiou| 


cx exp(-(d - Ax)'i'N~ r l (d - Ax) v /2) 


(6) 


199( 




Bouchet and Gispert 1998 


)) and the CMB power 



cx exp(- X 2 /2) 



(7) 



where N tt ' —< nn T > tt > is the noise covariance matrix. In 
this particular case, maximising V(x\y) amounts to find 



spectrum estimate. While potentially very interesting, this 
approach will not be generally applicable (at least effi- 
ciently) , and we now turn to a practical (general) solution 
of equation (||) by iterative means. 
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2.2. Practical implementation for large data sets 

We solve the map-making problem by adapting to our 
particular case the general "multi-grid method" (Press e1 
al. 1992| ) . Multi-grid methods are commonly used to speed 
up the convergence of a traditional relaxation method (in 
our case the Jacobi method, as in ( |Prunct et al. 2000| )) 
defined at resolution £ max (see below). A set of recursively 
defined coarser grids (£ < i m ax) are used as temporary 
computational space, in order to increase the convergence 
rate of the relaxation process. To be fully profitable, this 
algorithm implies for each resolution both a rebinning in 
space (resolution change) and in time (resampling). 

In this paper, we use the HEALPix pixelisation of the 
sphere (Gorski et al. 1998). In this scheme, the sphere 



is covered by 12 basic quadrilaterals, further divided re- 
cursively into pixels of equal area. The map resolution 
is labeled by N S ide- the number of pixels along the side 
of one basic quadrilateral. Hence, N S id e = 1 means that 
the sphere is covered by 12 large pixels only. The number 
of pixels is given by N pix = \2N 2 side . N side = 256 corre- 
sponds to a pixel size of 13.7 arcmin. For practical reasons, 
we need to define the resolution A: of a HEALPix map as 



N si de — 2 



(10) 



The "nested" pixel numbering scheme of HEALPix ( porski 
et al. 1998) allows an easy implementation of the coarsen- 
ing (k — > k — 1) and refining (k — ► k + 1) operators that 
we use intensively in our multigrid scheme. 

Let us now get into the details of our implementation 
and discuss successively the exact system we solve, the 
way we solve it and the actual steps of the multi-grid al- 
gorithm. 

2.2.1. Determining the working resolution k max 

We aim at solving for the optimal map &k at a given spatial 
resolution k using 



AlN~ l A k x^AlN- 1 d, 



(11) 



where A k is the "observation" operator (from spatial to 
temporal domain) and A\ is the "projection" operator 
(from temporal to spatial domain). In a noise- free exper- 
iment, the optimal map would be straightforwardly given 
by the co-added map (introducing the "co-addition" op- 
erator Pk) 



x k = P k d={AlA k )- l Ai d 



-i at 



(12) 



The time line is given by d — A^x^ + n where Xoo is the 
sky map at "infinite" resolution (fc = +oo in our nota- 
tions). In order to check the accuracy of this trivial noise- 
free map making, it is natural to compute the residual in 
the time domain with n = 



which will be non-zero in practice, as soon as one works 
with finite spatial resolution. We call this residual the pix- 
elisation noise. Since we assume here that the instrumen- 
tal beam is symmetric, the sky map is considered as the 
true sky convolved by, say, a Gaussian beam of angular 
diameter A8b- This introduces a low-pass spatial filter in 
the problem. In other words, as the resolution increase, the 
pixelisation noise should decrease towards zero. We have 
estimated that the order of magnitude of the pixelisation 
noise can be approximated by 



\\Pk\ 



Ag fc 



(14) 



The norms used in the above formula can be either the 
maximum over the time line (a very strong constraint) or 
the variance over the time line (a weaker constraint) . Since 
the pixelisation noise is strongly correlated with the sky 
signal, point sources or Galaxy crossings are potential can- 
didates for large and localised bursts of pixelisation noise. 
The correct working resolution k max is set by requiring 
that the pixelisation noise remains low compared to the 
actual instrumental noise, or equivalently 



A9 k . 



< 



A9 B 

S/N 



(15) 



Most of the CMB experiments are noise dominated along 
the time line, constraining the effective map resolution to 
be of the order of the instrumental beam or even larger. 
Note however that the pixelisation noise is strongly non 
Gaussian (point sources or Galaxy crossings) and can be 
always considered as a potential source of residual stripes 
in the final maps. 

2.2.2. The basic relaxation method: an approximate 
Jacobi solver 

Instead of solving for x we perform the change of variable 



y = x — P d 



(16) 



p k = A k x k 



(13) 



and solve instead for y which obeys 

PN- 1 Ay = PN- 1 {d-APd) or My = b (17) 

where we have multiplied each side of equation ( pd| ) by 
(A T A)~ 1 , the pixel hit counter. From now on, we also 
assume that the noise in the timestream is stationary 
and that its covariance matrix is normalised so that 
diag A^, 1 = /. The previous change of variable allows 
us to subtract the sky signal from the data: since we have 
chosen a resolution high enough to neglect the pixelisation 
noise, we have indeed APd ~ AqoXoo + APn and, conse- 
quently, d — APd ~ n — APn. The map making consists in 
two step: first compute a simple co-added map from the 
time line, and second, solve equation ( [i7| ) for the stripes 
map y. The final optimal map is obtained by adding these 
two maps. 



4 



Dore, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code 



It is worth mentioning that the stripes map is com- 
pletely independent of the sky signal, as soon as the pix- 
elisation noise can be ignored. It depends only on the scan- 
ning strategy through the matrix A and on the noise char- 
acteristics through the matrix N. Even if in principle this 
change of variable is irrelevant since it does not change the 
matrix to be inverted, it does in practice since d — APd 
is free from the fast variations of d (up to the pixelisation 
noise), as e.g. the Galaxy crossings, which are numerically 
damaging to the quality of the inversion. 

To solve for equation [l?], we follow the ideas of ( Prunet 
et al. I 000 ) and apply an approximate Jacobi relaxation 
scheme. The Jacobi method consists in the following re- 
laxation step to solve for our problem 



y 



" +1 = Ry n + D- l b and y° = 



(18) 



where D is the diagonal part of the matrix M and R is 
the residual matrix R = I — D^M . Computing the diag- 
onal elements of M — PN~ 1 A is rather prohibitive. The 
idea of Prunet et al. (2000) is to approximate D ~ I by 
neglecting the off-diagonal elements of iV _1 . The residual 
matrix then simplifies greatly and writes 



R = P(I - N~ 1 )A 



(19) 



The approximate Jacobi relaxation step is therefore de- 
fined as 



y 



n+l _ 



Ry n + b and y° = 



(20) 



One clearly sees that if this iterative scheme converges, 
it is towards the full solution of equation [l^. To perform 
these successive steps, it is extremely fruitful to remember 
the assumed stationarity of the noise. Indeed whereas this 
assumption implies a circulant noise covariance matrix in 
real space, it translates in Fourier space in the diagonality 
of the noise covariance matrix. This is naturally another 
formulation of the convolution theorem, since a stationary 
matrix acts on a vector as a convolution, and a diagonal 
matrix acts as a simple vector product, thus a convolution 
in real space is translated as a product in Fourier space. 
The point is that the manipulation of the matrix A^ 1 is 
considerably lighter and will be henceforth performed in 
Fourier space. Applying the matrix R to a map reduces 
then to the following operations in order 

1. "observe" the previous stripes map y n 

2. Fourier transform the resulting data stream 

3. apply the low-pass filter W = I — N^ 1 

4. inverse Fourier transform back to real space 

5. co-add the resulting data stream into the map y n+1 . 

Assuming that the normalised noise power spectrum 
can be approximated by P(f) = 1 + (/o//)°, the low-pass 
filter associated to each relaxation step is given by 



W(J) = 



fa 
JO 



Since both A and P are norm-preserving operators, the 
norm of the increment Ay n — y n — y n_1 between step 
nandn+1 decreases as ||Ay™ +1 || < W(f mm )\\Ay n \\, 
where f m in is a minimal frequency in the problem. Since 
W(fmin) < 1, we see that the approximate Jacobi relax- 
ation scheme will converge in every case, which is good 
news. On the other hand, since W(f m i n ) ~ 1, the ac- 
tual convergence rate of the scheme is likely to be very, 
very slow, which is bad news (cf. figure ^ for a graphi- 
cal illustration) . The fact that this algorithm is robust, 
but dramatically slow is a well-known property of the Ja- 
cobi method. The multi-grid method is also well known 
to solve this convergence speed problem. Note that if the 
convergence is reached, the solution we get is the optimal 
solution, i.e. similar to the one that would be obtained by 
a full matrix inversion. 



2.2.3. Multigrid relaxation 

The multi-grid method for solving linear problems is de- 



scribed in great details in (Press et al. 1992). At each 
relaxation step at level k — k max , our target resolution, 
we can define the error ejj = Vk ~ Vk and the residual 
r£ = Mkj)% — bk- Both are related through 



M k e n k = rt 



(22) 



IS + f a 



(21) 



If we are able to solve exactly for equation (|22j) , the overall 
problem is solved. The idea of the multi-grid algorithm is 
to solve approximately for equation (|2^) using a coarser 
"grid" at resolution k — 1, where the relaxation scheme 
should converge faster. We thus need to define a fine-to- 
coarse operator, in order to define the new problem on 
the coarse grid and solve for it. We also need a coarse-to- 
fine operator in order to inject the solution onto the fine 
grid. The approximation to the error e£ is finally added 
to the solution. The coarse grid solver is usually applied 
after a few iterations on the fine level have been performed 
(in practice we perform 3 to 5 iterations). Naturally, since 
the solution to the problem on the coarse level relies also 
on the same relaxation scheme, it can be itself acceler- 
ated using an even coarser grid. This naturally leads to a 
recursive approach of the problem. 

We defined our fine-to-coarse operator to be an aver- 
aging of the values of the 4 fine pixels contained in each 
coarse pixel. The coarse-to-fine operator is just a straight 
injection of the value of the coarse pixel into the 4 fine 
pixels. The most important issue is the temporal rebin- 
ning of the data stream since the speed of the iterative 
scheme at a given level is actually set by the two Fourier 
transforms. We performed that resampling by simply tak- 
ing each time we go up one level one point out of two. At 
the lower resolutions, this reduction is such that the itera- 
tion cost is negligible when compared to that at higher k; 
it allows fast enough iterations that full convergence may 
be reached. In practice we choose a minimal level k = 3 
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defined by N S id e = 8 and iterate a few hundred times to 
reach exact convergence. 



Finally the navigation through levels allows several op- 



We introduced 5 distinct levels of resolution defined by 
their N S id e parameter in the HEALPix package ( porskj| 
et al. 1998). The higher resolution level is imposed by 



tions to be taken. Either we go up and down through all 
the levels successively (the so-called "V-cycle" ) or we fol- 
low more intricate paths {e.g. the "W-cycle" where at a 
given level we go down and up all the lower levels before 
going up and vice- versa). Since it turns out that levels are 
relatively disconnected in the sense that the scales well 
solved at a given resolution are only slightly affected by 
the solution at a higher level, the "V-cycle" is appropriate 
and is the solution we adopt in the two following configu- 
rations. 

3. Practical application to ARCHEOPS and 
TopHat experiments 

We now aim at demonstrating the capabilities of this al- 
gorithm with simulated data whose characteristics are in- 
spired by those of the ARCHEOPS and TopHat experi- 
ments. 

3.1. Simulation 

The ARCHEOPS and TopHat experiments are very sim- 
ilar with respect to their scanning strategy. Indeed both 
use a telescope that simply spins at a constant rate (re- 
spectively 3 and 4 RPM) about its vertical axis. Thus due 
to Earth rotation the sky coverage of both is performed 
through scan circles whose axis slowly drifts on the sky. 
Nevertheless, because of the different launch points (re- 
spectively on the Arctic circle (Kiruna, Sweden) and in 
Antarctica (McMurdo)) and their different telescope axis 
elevation (respectively ~ 41° or 12°) they do not have the 
same sky coverage. 

Otherwise the two experiments do not use the same 
bolometers technology, neither the same bands nor have 
the same number of bolometers. But even if we try to be 
fairly realistic, our goal though is not to compare their re- 
spective performances but rather to demonstrate two ap- 
plications of our technique in different settings. We then 
simulate for each a data stream of ~ 24 hrs duration 
with respectively a sampling frequency of 171 and 64 Hz. 
The TODs contain realistic CMB and Galactic signal for 
a band at 143 GHz. Note that this is a one day data 
stream for TopHat (out of 10 expected) and that this fre- 
quency choice is more appropriate for ARCHEOPS than 
for TopHat (whose equivalent band is around 156 GHz), 
but this is realistic enough for our purpose. We generated 
a Gaussian noise time stream with the following power 
spectrum P(f) oc (1 + (f knee / f)")- 1 with f knee = 0.24 
and 1 Hz, and a — 1.68 and 1. The noise amplitude per 
channel is choosen so that it corresponds for ARCHEOPS 
(24 hours) and TopHatQ (10 days of flight) to 30/8 fiK on 
average per 20' pixel, with a beam FWHM of 10' / 20'. 

3 Lloyd Knox private communication 



the pixclisation noise level requirement (section 2.2.3| ) to 
N S ide = 256 (pixel size ~ 13.7') whereas the lower one is 
Nside — 8 (pixel size ~ 7.3°). Therefore these two config- 
urations each offer an interesting test since they differ by 
the sky coverage and the noise power spectrum. We iter- 
ate 3 times at each level except at the lowest one where 
we iterate 100 times. 



3.2. Results 

The algorithm is as efficient in both situations. Whereas 
for ARCHEOPS whose timeline is longer due to the higher 
sampling frequency it took 2.25 hours on a SGI ORIGIN 
2000 single processor, it took around 1.5 hours for the 
TopHat daily data stream. 

In figures [j] and || we depict from top to bottom and 
from left to right the initial co-added data, the recon- 
structed noise map, the hit map, i.e. the number of hit 
per pixels at the highest worked out resolution, the initial 
signal map as well as the reconstructed signal map and the 
error map. We see that the destriping is excellent in both 
situations and the signal maps recovered only contain an 
apparently isotropic noise. We note the natural correla- 
tion between the error map and the hit map. Finally, we 
stress that no previous filtering at all has been applied. 

Figure |^ shows how the noise map is reconstructed at 
various scales. This is a mere illustration of our multi-grid 
work. 



3.3. Tests 

At this point, some tests are required to cautiously vali- 
date our method. First, as was stated below, as soon as 
the iterative algorithm has converged, the solution is by 
construction the optimal solution, similar to the one that 
would be obtained by the full matrix inversion. As a cri- 
terium for convergence we required the 2-norm of residuals 
to be of the order of the machine accuracy. 

We initially have a Gaussian random noise stream fully 
characterised by its power spectrum. Therefore an impor- 
tant test is to check wether the deduced noise stream (by 
"observing" the stripe map, see § 4.3 for further details 



definition) is Gaussian and has the same power spectrum. 
On figure^ we ensure that the evaluated noise time stream 
is indeed gaussian. As depicted on figure ||, where we plot 
in the ARCHEOPS case both the analytical expression of 
the spectrum according to which we generate the timeline 
and its recovered form, the agreement is excellent. We re- 
call that we assumed at the beginning a perfect knowledge 
of the noise in order to define the filters. This is naturally 
unreali stic but the issue of noise evaluation is discussed in 
section 4J3 below. We plotted as well the probability dis- 
tribution function (PDF) of the final error map, i.e. the 
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Fig. 1. Simulated ARCHEOPS Kiruna flight : From top to bottom and from left to right, the co-added map and the 
input Galactic + CMB signal map, the reconstructed noise (stripes) and signal maps, the hit count and the error map 
(Arbitrary unit). The fact that the coverage is not fully symmetrical is due to the fact that we considered slightly less 
than 24hr. Mollweide projection with pixels of 13.7' (HEALPix N s id e = 256). Arbitrary units. 



Fig. 2. Simulated TopHat one day flight : From top to bottom and from left to right, the co- added map and the 
input CMB + Galaxy signal map, the reconstructed noise (or stripes) and signal map, the hit count and error map. 
The fact that the coverage is not fully symmetrical is due to the fact that we consider only 18.2hr of flight. Gnomonic 
projection with pixel of 13.7' (HEALPix N S id e — 256). Note that the slight visible stripping is correlated to the 
incomplete rotation pattern. Arbitrary units. 



Fig. 3. Multi-grid noise map recovery : In this plot we show in the ARCHEOPS Kiruna case how the noise map is 
reconstructed at various levels, corresponding respectively to N S id e — 8, 16, 32, 64. 



recovered noisy signal map minus the input signal map 
(figure |j). This PDF is well fitted by a gaussian whose 
parameters are given in figure f|. The PDF of the error 
map displays some non-gaussian wings. Let us recall here 
that this is no surprise here because of the non-uniform 
sky coverage as well as the slight residual stripping, both 
due to the non-ideal scanning strategy, i.e. that produces 
a non-uniform white noise in pixel space. 

Another particularly important test consists in check- 
ing the absence of correlation between the recovered noise 
map and the initial signal map. We could not find any 
which is no surprise since we are iterating on a noise map 
(see § 2.2.2 ) which does not contain any signal up to the 



pixelisation noise, that is ensured to be negligible with re- 
gards to the instrumental noise by choice of the resolution 



kn 



(see § 2.2.1) 



4. Discussion 

4-.1. Why is such an algorithm so efficient ? 

The efficiency of such a method can be intuitively un- 
derstood. Indeed, although the Jacobi method is known 
to converge very safely it suffers intrinsically from a very 
slow convergence for large-scale correlations (which origi- 
nate mainly in the off diagonal terms of the matrix to be 
inverted) (Press et al. 1992). This is illustrated on figure 
|6|: there we show the maps of residuals after solving this 
system using a standard Jacobi method on simulated data 
with 50, 100, 150, and 200 iterations. We used the same 
simulation and therefore the same sky coverage. Obviously 
the largest structures are the slowest to converge (imply- 
ing observed large scale residual patterns). As a conse- 
quence it seems very natural to help the convergence by 
solving the problem at larger scales. Whereas large-scale 
structures will not be affected by a scale change, smaller 
structures will converge faster. 



4-2. Scaling s 

We have found that this multi-grid algorithm translates 
naturally in a speed up greater than 10 as compared to a 
standard Jacobi method. This is illustrated in figure (J7|) 
where we plotted the evolution of the 2-norm of residu- 
als for the two methods in terms of the number of iter- 
ations in 'cycle units'. One cycle corresponds to 8 itera- 
tions at level k max for a standard Jacobi method whereas 
it incorporates addionally going up and down all the low- 
est levels in the multi-grid case. Thus the cycle timing 
is not exactly the same but the difference is negligible 
since the limiting stages are definitely the iterations per- 
formed at maximum resolution. Note the fact that the effi- 
ciency of the multi-grid method allows us to solve exactly 
the system up to the machine accuracy (small rebounds 
at the bottom of the curve) in approximately 135mn for 
{Ntod ~ 8 10 6 ,AAp« - 8 10 5 ) on a SGI ORIGIN 2000 sin- 
gle processor. Since the limiting stages are the FFT's at 
the higher levels, this algorithm scales as 0(J\f to d hi Ntod)- 
In terms of memory storage it scales naturally as 0{Mtod) 
since one crucial feature of this iterative method is to han- 
dle only vectors and never an explicit matrix. The scaling 
of this problem is formally independent of the number of 
pixels Afpi x . 

4-3. The noise estimation issue and noise covariance 
matrix estimation 

The estimation of the statistical properties of the noise in 
the timestream is an important issue for this kind of algo- 
rithm. Indeed, whereas till now we have assumed a perfect 
prior knowledge of the noise power spectrum in order to 
define the filters, it might not be that easy in practice since 
we can not separate the signal from the noise. We will 
therefore aim at making a joint estimation of the signal 
|and the noise. This has been pioneered recently by (Fer 



Dore, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code 



7 



Fig. 6. Residual map after 50, 100, 150 and 200 iterations of a standard Jacobi method. This has been performed on 
simulated data for one bolometer with a nominal noise level. The sky coverage is that of ARCHEOPS coming Kiruna 
flight. The residual large scale patterns illustrate the difficulties the standard Jacobi method faces to converge. The 
stripes free area are the ones of scan crossing (see the hit map in figure 111) . 



PDF of evaluated noise timestrean- 1 




Frequency [Hz] 



PDF of error mop 




Fig. 4. In the TopHat case we plot from top to bottom 
the recovered probability distribution function of the noise 
stream evaluated along the timeline as well as the error 
map PDF. In this two cases a fit to a gaussian has been 
performed whose parameters are written inside the figures. 
No significant departure from gaussianity are detected. 
The arbitrary units are the same as the ones used for the 
previously shown maps. 



reira and Jaffc 2000) and implemented independently by 
(Prunet et al. 2000). The latter implementation is rather 
straightforward in our particular case since it just implies 
reevaluating the filters after a number of iterations, given 
the (current) estimation of the signal map and thus of 
the signal data stream. Nevertheless its non-obvious con- 
vergence properties have to be studied carefully through 



Fig. 5. Recovery of the noise power spectrum in the 
Archeops case (top: linear x-axis, bottom: log x-axis): The 
red thin dashed line shows the initial analytic noise power 
spectrum used to generate the input noise stream and the 
blue thick line denotes the recovered one after 6 iterations. 
The recovered one has been binned as described in section 
4.3 and both have been normalised so that the white high 
frequency noise spectrum is equal to one. The agreement 
is obviously very good. No apparent bias is visible. Note 
that a perfect noise power spectrum prior knowledge has 
been used in this application. 



simulations. Making use of ([!(]) our evaluation of the noise 



timeline fi n at the 



-,th 



iteration and at level max is 



n n =d-A{y k n ma +Pd) 



(23) 
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Fig. 7. The evolution of the 2-norm of residuals with 
the number of iterations at level max. Whereas the blue 
dashed line is standard iterative Jacobi, the solid red line 
is the iterative multi-grid Jacobi method. A full multi- 
grid cycle incorporates 3 iterations at level max before 
going down and up all the levels. The sharp jumps corre- 
sponds to the moment when the scheme reach again level 
max and thus take the full benefits from the resolution at 
lower levels. Note that the very sharp features after ~ 100 
iterations are due to the fact we reached the machine accu- 
racy which is almost out of reach for a standard iterative 
method. 



Then we compute its spectrum and (re-) define the re- 
quired filters. We then go through several multi-grid 
cycles (5 in the above demonstrated case) before re- 
evaluating the noise stream. Very few evaluations of the 
noise are needed before getting a converged power spec- 
trum (around 2). In such an implementation, no noise 
priors at all are assumed. This is illustrated on one par- 
ticular worked out example in the case of a 4 hours 
ARCHEOPS like flight (more detailed considerations will 
be discussed somewhere else). To reduce the number of 
degrees of freedom we bin the evaluated noise power spec- 
trum using a constant logarithmic binning (A In / = 0.15 
in our case) for / < 2 fknee and a constant linear bin- 
ning (A / = 0.08 Hz in our case) for higher frequency. 
The figure || shows the genuine and evaluated noise power 
spectrum. The initial noise power spectrum was a real- 
istic one P(f) oc (1 + (fknee/ '/)") to which we added 
some small perturbations (the two visible bumps) to test 
the method. Note the small bias around the telescope spin 
frequency at f sp in = 0.05Hz: this is illustrative of the diffi- 
culties we fundamentally face to separate signal and noise 
at this particular frequency through equation ( |l6| ) . Natu- 
rally, this bias was not present in the case demonstrated on 
figure U where we assumed a prior knowledge of the spec- 
trum. This possible bias forced us to work with a coarser 
binning (A In/ = 1.) in the 1// part of the spectrum 



till the convergence is reached, i.e. we evaluate the noise 
power spectrum with the previously mentioned binning 
only at the last step. Proceeding this way, the conver- 
gence towards the correct spectrum is both fast (3 noise 
evaluations) and stable. 

Second, the output of any map-making should contain 
as well an evaluation of the map noise covariance matrix 
(A T N~ 1 A)~ 1 . Given such a fast algorithm and given an 
evaluation of the power spectrum, it is natural to obtain 
it through a Monte-Carlo algorithm fueled with various 
realizations of the noise timeline generated using the eval- 
uated power spectrum. This part will be presented in a 
future work. However we illustrate it very briefly by a 
very rough Ci determination (which is in no way an ap- 
propriate Cg estimate) . To this purpose we perform a one 
day TopHat like simulation including only the CMB sig- 
nal plus the noise. From this data stream we obtain an 
"optimal" signal map as well as an evaluation of the noise 
power spectrum using the previously described algorithm. 
With the help of the anaf ast routine of the HEALPix 
package we calculate this way a rough (jf' 9nal . Using the 
estimated noise power spectrum we generate 10 realisa- 
tions of the noise and get consequently 10 "optimal" noise 
maps. For each of them we measure as before Cg and aver- 
age them to obtain C™ mse . In order to debiase the signal 
power spectrum recovered in this way, we substract C™ mse 
to (J^ 9nal . The power spectrum obtained in this manner 
includes as well some spatial correlations due to some light 
residual stripping and thus does not correspond to white 
noise (at least in the low I part) . For the aim of compar- 
ison we compute the power spectrum of the input signal 
measured the same way, C l ™ pwt , and plotted both C™ put 
as well as Cg l9na — C r g mse averaged in linear bands of 
constant width A£ = 80. The agreement is obviously very 
good as illustrated on figure ||| The error bars take into 
account both the sampling in duced variance as well as the 
beam spreading ( Knox 1995| ) . A few comments need to be 
made at this point. This is in no way an appropriate Cg 
measurement since we are not dealing properly here with 
the sky coverage induced window function which triggers 
some spurious correlations within the Cgs. We thus do 
not take into account the full-structure of the map noise 
covariance matrix. Nevertheless, the efficiency of such a 
rough method is encouraging for more detailed implemen- 
tation and a full handling of the noise covariance matrix. 

Note finally that this is a naturally parallelised task 
which should therefore be feasible very quickly. 

4-4- Application to genuine data and hypothesis 

The application to genuine data could be problematic if 
our key hypothesis were not to be fulfilled thus we have 
to discuss them. Concerning the noise, we assumed that 
it is Gaussian in order to derive our map estimator and 
stationary in order to exploit the diagonality of its noise 



Dore, Teyssier, Bouchet, Vibert & Prunet: MAPCUMBA: a fast multi-grid CMB map-making code 



9 




Fig. 9. In the case of a one day flight of the TopHat ex- 
periment we perform a very approximate evaluation of 
the recovered signal band powers (black) which has to be 
compared to the input signal power spectrum (red line). 
Both have been measured the same way using the anaf ast 
routine of the HEALPix package. These band powers 
(A£ = 80) have been performed using a fast Monte-Carlo 
evaluation of C™ mse (dashed blue line) which does not cor- 
respond exactly to white noise since there remains some 
spatial correlations. This constitutes in no way an appro- 
priate Ci measurement but is an encouraging step towards 
a full Monte-Carlo treatment. 



l.E-01 I 

0.01 



Fig. 8. Evaluation of the noise power spectrum (top: lin- 
ear x-axis, bottom: log x-axis): the red thin dashed line 
shows the initial noise power spectrum obtained from the 
input noise stream and the blue thick line denotes the re- 
covered one after 5 iterations. Both have been smoothed 
and normalised so that the white high frequency noise 
spectrum is equal to one. The agreement is obviously very 
good. Note that no noise priors at all have been used in 
this evaluation. 



covariance matrix in Fourier space. Both hypothesis are 
reasonable for a nominal instrumental noise, at least par- 
tially on (sufficiently long) pieces of timeline. If not, we 
would have to cut the bad parts and replace them by a 
constrained realization of the noise in these parts. Con- 
cerning the signal, no particular assumptions are needed 
in the method we are presenting. At this level, we ne- 
glected as well the effect of the non perfect symmetry of 
the instrumental beam. This effect should be quantified 

Pore et al. 2000| ). 



for a real experiment (Wu et al. 2000 



Another technical hypothesis is the negligibility of the pix- 
elisation noise with respect to the instrumental noise but 
since this is fully under control of the user it should not 
be a problem. 



In this paper, we have presented an original fast map- 
making method based on a multi-grid Jacobi algorithm. 
It naturally entails the possibility of a joint noise/signal 
estimation. We propose a way to generate the noise co- 
variance matrix and illustrate its ability on a simple Ce 
estimation. The efficiency of this method has been demon- 
strated in the case of two coming balloon experiments, 
namely ARCHEOPS and TopHat but it naturally has a 
more general range of application. This tool should be of 
natural interest for Planck and this will be demonstrated 
somewhere else. Furthermore, due to the analogy of the 
formalisms, this should have some applications as well in 
the component separation problem. 

We hope to apply it very soon on 
ARCHEOPS data. The FORTRAN 90 code whose 
results have been presented is available at 
http : //ulysse . iap . f r/download/mapcumba/. 
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